Doubly Truncated Weibull Minimum Distribution (truncweibull_min)#
scipy.stats.truncweibull_min is the Weibull minimum distribution restricted to an interval and re-normalized.
Truncation is common when you only observe values within a window (instrument limits, study design, reporting rules) but the underlying phenomenon is still reasonably Weibull.
Learning goals#
Write down the PDF/CDF/PPF and understand how truncation re-normalizes a base Weibull.
Compute moments using incomplete gamma functions and verify them numerically.
Interpret how the parameters
(c, a, b, loc, scale)change the shape and the support.Sample from the distribution with a NumPy-only inverse-transform algorithm.
Use
scipy.stats.truncweibull_minfor evaluation, simulation, and fitting.
Notebook roadmap#
Title & classification
Intuition & motivation
Formal definition (PDF/CDF/PPF)
Moments & properties
Parameter interpretation
Derivations (expectation, variance, likelihood)
Sampling & simulation (NumPy-only)
Visualization (PDF, CDF, Monte Carlo)
SciPy integration (
scipy.stats.truncweibull_min)Statistical use cases
Pitfalls
Summary
Prerequisites#
comfort with PDF/CDF and expectation
basic calculus substitutions
familiarity with the Gamma / incomplete Gamma functions (helpful but not required)
import numpy as np
import plotly
import plotly.express as px
import plotly.graph_objects as go
import os
import plotly.io as pio
from plotly.subplots import make_subplots
import scipy
from scipy import integrate, optimize, special, stats
from scipy.stats import truncweibull_min
pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
np.set_printoptions(precision=6, suppress=True)
rng = np.random.default_rng(42)
# Record versions for reproducibility.
VERSIONS = {"numpy": np.__version__, "scipy": scipy.__version__, "plotly": plotly.__version__}
VERSIONS
{'numpy': '1.26.2', 'scipy': '1.15.0', 'plotly': '6.5.2'}
1) Title & Classification#
Item |
Value |
|---|---|
Name |
|
Type |
Continuous |
Base distribution |
Weibull minimum ( |
Support (standardized) |
\(x \in (a, b]\) |
Parameter space |
\(c>0\), \(0\le a<b\le\infty\) |
Location/scale |
|
SciPy defines the truncation bounds in standardized form:
so the support in the original units is
Unless stated otherwise, we work in the standardized form (loc=0, scale=1).
2) Intuition & Motivation#
What it models#
The Weibull minimum distribution is a staple model for time-to-failure / lifetime phenomena. A truncated Weibull is appropriate when the physical process is still Weibull-like, but you only observe outcomes inside a window:
Left truncation (e.g., detection limits): values below a threshold are unobservable.
Right truncation (e.g., instrument saturation, administrative cutoffs): values above a threshold are not recorded.
Double truncation: both happen.
Truncation is not the same as censoring: with truncation, observations outside the window never appear in the dataset.
Typical real-world use cases#
Reliability testing with a fixed observation period (components only recorded if they fail between times \(u_\ell\) and \(u_r\)).
Survival analysis with delayed entry (left truncation) and administrative endpoints (right truncation).
Environmental measurements (wind speeds, rainfall intensities) where sensors have detection/saturation limits.
Relations to other distributions#
As \(a\to 0\) and \(b\to\infty\),
truncweibull_minreduces to the standard Weibull minimum distribution.When \(c=1\), Weibull becomes Exponential, and
truncweibull_minbecomes a truncated exponential.When \(c=2\) (with a suitable scale), Weibull relates to the Rayleigh distribution.
It is a special case of a generic truncated distribution: conditioning a base distribution to lie in \((a,b]\).
3) Formal Definition#
Base Weibull (standardized)#
For shape parameter \(c>0\), the (standardized) Weibull minimum distribution has
Its survival function is \(S_W(x;c)=e^{-x^c}\).
Truncated Weibull minimum#
Define the normalization constant
Then the truncated density is
The CDF is
A convenient form for the quantile function (PPF) is obtained by solving for \(e^{-x^c}\):
Location/scale#
If \(Y\sim\mathrm{truncweibull\_min}(c,a,b)\) in standardized form, then
has support \((\mathrm{loc}+a\,\mathrm{scale},\; \mathrm{loc}+b\,\mathrm{scale}]\) and
def _logZ_truncweibull_min(c: float, a: float, b: float) -> float:
"""Stable log normalizer log(exp(-a^c) - exp(-b^c)).
Notes
-----
For b = inf, this is simply -a^c.
"""
if not (c > 0):
raise ValueError("c must be > 0")
if not (a >= 0):
raise ValueError("a must be >= 0")
if not (b > a):
raise ValueError("b must be > a")
A = a**c
if np.isinf(b):
return -A
B = b**c
# exp(-A) - exp(-B) = exp(-A) * (1 - exp(-(B-A)))
return -A + np.log1p(-np.exp(-(B - A)))
def truncweibull_min_logpdf(x: np.ndarray, c: float, a: float, b: float) -> np.ndarray:
"""Log-PDF of the standardized truncweibull_min(c, a, b)."""
x = np.asarray(x, dtype=float)
logpdf = np.full_like(x, fill_value=-np.inf)
mask = (x > a) & (x <= b)
if not np.any(mask):
return logpdf
logZ = _logZ_truncweibull_min(c, a, b)
xm = x[mask]
logpdf[mask] = np.log(c) + (c - 1.0) * np.log(xm) - xm**c - logZ
return logpdf
def truncweibull_min_pdf(x: np.ndarray, c: float, a: float, b: float) -> np.ndarray:
"""PDF of the standardized truncweibull_min(c, a, b)."""
return np.exp(truncweibull_min_logpdf(x, c, a, b))
def truncweibull_min_cdf(x: np.ndarray, c: float, a: float, b: float) -> np.ndarray:
"""CDF of the standardized truncweibull_min(c, a, b)."""
x = np.asarray(x, dtype=float)
out = np.zeros_like(x)
out[x > b] = 1.0
mask = (x > a) & (x <= b)
if not np.any(mask):
return out
A = a**c
xm = x[mask]
delta_x = xm**c - A
# Numerator: exp(-a^c) - exp(-x^c) = exp(-a^c) * (1 - exp(-(x^c-a^c)))
num = -np.expm1(-delta_x)
if np.isinf(b):
out[mask] = num
return out
delta_b = b**c - A
den = -np.expm1(-delta_b)
out[mask] = num / den
return out
def truncweibull_min_ppf(q: np.ndarray, c: float, a: float, b: float) -> np.ndarray:
"""PPF (quantile function) of the standardized truncweibull_min(c, a, b)."""
q = np.asarray(q, dtype=float)
if np.any((q < 0) | (q > 1)):
raise ValueError("q must be in [0, 1]")
A = a**c
sa = np.exp(-A)
sb = 0.0 if np.isinf(b) else np.exp(-(b**c))
s = (1.0 - q) * sa + q * sb
return (-np.log(s)) ** (1.0 / c)
# Sanity check against SciPy
c0, a0, b0 = 2.5, 0.3, 2.0
x = np.linspace(a0 + 1e-6, b0, 500)
q = np.linspace(1e-6, 1 - 1e-6, 500)
pdf_err = np.max(np.abs(truncweibull_min_pdf(x, c0, a0, b0) - truncweibull_min.pdf(x, c0, a0, b0)))
cdf_err = np.max(np.abs(truncweibull_min_cdf(x, c0, a0, b0) - truncweibull_min.cdf(x, c0, a0, b0)))
ppf_err = np.max(np.abs(truncweibull_min_ppf(q, c0, a0, b0) - truncweibull_min.ppf(q, c0, a0, b0)))
pdf_err, cdf_err, ppf_err
(3.3306690738754696e-16, 2.220446049250313e-16, 0.0)
4) Moments & Properties#
Raw moments#
For \(k\ge 0\), the \(k\)-th raw moment exists (and is finite) whenever the support avoids \(0\) or the left tail is integrable. For the standardized distribution with truncation \((a,b]\) and \(a\ge 0\):
With the substitution \(t=x^c\), this becomes an incomplete gamma integral:
where \(\gamma(s,x)\) is the lower incomplete gamma function.
Mean / variance / skewness / kurtosis#
Let \(m_k = \mathbb{E}[X^k]\). Then
Mean: \(\mu=m_1\)
Variance: \(\sigma^2 = m_2 - m_1^2\)
Skewness: \(\gamma_1 = \mathbb{E}[(X-\mu)^3]/\sigma^3\)
Excess kurtosis: \(\gamma_2 = \mathbb{E}[(X-\mu)^4]/\sigma^4 - 3\)
MGF / characteristic function#
The moment generating function (MGF) is
If \(b<\infty\), then the support is bounded and the MGF exists for all real \(t\).
If \(b=\infty\), the MGF existence matches the base Weibull (e.g. for \(c=1\) it only exists for \(t<1\)).
The characteristic function always exists:
Closed forms are generally not simple; numerical quadrature is typical.
Entropy#
The differential entropy is
SciPy provides truncweibull_min.entropy, and we can also compute it numerically.
def truncweibull_min_raw_moment(k: float, c: float, a: float, b: float) -> float:
"""Raw moment E[X^k] in standardized form.
Uses the formula with lower incomplete gamma:
E[X^k] = (γ(1+k/c, b^c) - γ(1+k/c, a^c)) / (exp(-a^c) - exp(-b^c)).
"""
if k < 0:
raise ValueError("This helper is written for k>=0.")
s = 1.0 + k / c
A = a**c
B = np.inf if np.isinf(b) else b**c
gamma_s = special.gamma(s)
num = gamma_s * (special.gammainc(s, B) - special.gammainc(s, A))
den = np.exp(-A) - (0.0 if np.isinf(B) else np.exp(-B))
return float(num / den)
def truncweibull_min_moments(c: float, a: float, b: float) -> dict:
m1 = truncweibull_min_raw_moment(1, c, a, b)
m2 = truncweibull_min_raw_moment(2, c, a, b)
m3 = truncweibull_min_raw_moment(3, c, a, b)
m4 = truncweibull_min_raw_moment(4, c, a, b)
var = m2 - m1**2
mu3 = m3 - 3 * m1 * m2 + 2 * m1**3
mu4 = m4 - 4 * m1 * m3 + 6 * (m1**2) * m2 - 3 * m1**4
skew = mu3 / (var ** 1.5)
kurt_excess = mu4 / (var**2) - 3.0
return {"mean": m1, "var": var, "skew": skew, "kurt_excess": kurt_excess}
c1, a1, b1 = 1.7, 0.2, 2.0
ours = truncweibull_min_moments(c1, a1, b1)
scipy_mean, scipy_var, scipy_skew, scipy_kurt = truncweibull_min.stats(c1, a1, b1, moments="mvsk")
ours, (scipy_mean, scipy_var, scipy_skew, scipy_kurt)
({'mean': 0.8838055783319058,
'var': 0.19039626228275253,
'skew': 0.5000215465747762,
'kurt_excess': -0.585680805798487},
(0.8838055783319058,
0.19039626228275253,
0.5000215465747815,
-0.585680805798499))
def mgf_numeric(t: float, c: float, a: float, b: float) -> float:
"""Numerical MGF M(t)=E[e^{tX}] via quadrature."""
def integrand(x: float) -> float:
return float(np.exp(t * x) * np.exp(truncweibull_min_logpdf(np.array([x]), c, a, b)[0]))
val, err = integrate.quad(integrand, a, b, limit=400)
return float(val)
def cf_numeric(t: float, c: float, a: float, b: float) -> complex:
"""Numerical characteristic function φ(t)=E[e^{itX}] via quadrature."""
def integrand_cos(x: float) -> float:
lp = truncweibull_min_logpdf(np.array([x]), c, a, b)[0]
return float(np.cos(t * x) * np.exp(lp))
def integrand_sin(x: float) -> float:
lp = truncweibull_min_logpdf(np.array([x]), c, a, b)[0]
return float(np.sin(t * x) * np.exp(lp))
re, _ = integrate.quad(integrand_cos, a, b, limit=400)
im, _ = integrate.quad(integrand_sin, a, b, limit=400)
return complex(re, im)
def entropy_numeric(c: float, a: float, b: float) -> float:
"""Numerical differential entropy -E[log f(X)]."""
def integrand(x: float) -> float:
lp = truncweibull_min_logpdf(np.array([x]), c, a, b)[0]
if not np.isfinite(lp):
return 0.0
return float(-np.exp(lp) * lp)
val, _ = integrate.quad(integrand, a, b, limit=400)
return float(val)
t_values = [-2.0, -1.0, 0.5]
mgf_vals = [mgf_numeric(t, c1, a1, b1) for t in t_values]
cf_vals = [cf_numeric(t, c1, a1, b1) for t in t_values]
H_scipy = float(truncweibull_min.entropy(c1, a1, b1))
H_num = entropy_numeric(c1, a1, b1)
mgf_vals, cf_vals, H_scipy, H_num
([0.23492536535757963, 0.45106303854527163, 1.5944101662659167],
[(-0.08811862717309919-0.6682569212463516j),
(0.5811634487536085-0.698157117953174j),
(0.8829961202882072+0.4168074733874823j)],
0.4714474041561004,
0.4714474041561004)
5) Parameter Interpretation#
Shape parameter c (Weibull shape)#
In the (untruncated) Weibull minimum distribution:
\(c<1\) produces a decreasing density on \((0,\infty)\) with a heavy right tail (sub-exponential).
\(c=1\) reduces to the exponential distribution.
\(c>1\) produces a density that rises from 0 and then decays; the hazard rate is increasing.
After truncation, the same tendencies remain, but the distribution is confined to \((a,b]\).
Truncation parameters a and b#
Increasing
aremoves the left part of the support and shifts mass to the right.Decreasing
bremoves the upper tail and forces the distribution to concentrate belowb.Limits:
\(b\to\infty\) gives a left-truncated Weibull.
\(a\to 0\) and \(b\to\infty\) recovers the base Weibull.
loc and scale#
locshifts the distribution.scalestretches/compresses it.The truncation bounds in the original units are \(u_\ell = \mathrm{loc}+a\,\mathrm{scale}\) and \(u_r = \mathrm{loc}+b\,\mathrm{scale}\).
def plot_pdf_family_by_c(a: float, b: float, c_values: list[float]) -> go.Figure:
x = np.linspace(a + 1e-6, b, 800)
fig = go.Figure()
for c in c_values:
fig.add_trace(go.Scatter(x=x, y=truncweibull_min_pdf(x, c, a, b), mode="lines", name=f"c={c}"))
fig.update_layout(
title="PDF changes with shape c (fixed a,b)",
xaxis_title="x",
yaxis_title="density",
legend_title="shape",
)
return fig
def plot_pdf_family_by_trunc(c: float, trunc_pairs: list[tuple[float, float]]) -> go.Figure:
fig = go.Figure()
for a, b in trunc_pairs:
x = np.linspace(a + 1e-6, b, 800)
fig.add_trace(go.Scatter(x=x, y=truncweibull_min_pdf(x, c, a, b), mode="lines", name=f"a={a}, b={b}"))
fig.update_layout(
title="PDF changes with truncation (fixed c)",
xaxis_title="x",
yaxis_title="density",
legend_title="(a,b)",
)
return fig
fig1 = plot_pdf_family_by_c(a=0.2, b=2.0, c_values=[0.7, 1.0, 1.7, 3.0])
fig2 = plot_pdf_family_by_trunc(c=1.7, trunc_pairs=[(0.0, 2.0), (0.2, 2.0), (0.2, 1.2), (0.8, 2.0)])
fig1.show()
fig2.show()
6) Derivations#
(a) Raw moment and expectation#
Let \(X\sim\mathrm{truncweibull\_min}(c,a,b)\) in standardized form. For \(k\ge 0\):
Substitute \(t=x^c\) so that \(dt=c x^{c-1}dx\) and \(x^k=t^{k/c}\):
Dividing by \(Z(c,a,b)=e^{-a^c}-e^{-b^c}\) yields the raw moment formula used above.
(b) Variance#
Once \(m_1=\mathbb{E}[X]\) and \(m_2=\mathbb{E}[X^2]\) are available,
(c) Likelihood#
Given i.i.d. observations \(x_1,\dots,x_n \in (a,b]\), the log-likelihood (standardized form) is
With loc/scale, use \(y_i=(x_i-\mathrm{loc})/\mathrm{scale}\) and add the Jacobian term \(-n\log(\mathrm{scale})\).
def loglike_c_given_ab(c: float, x: np.ndarray, a: float, b: float) -> float:
"""Log-likelihood as a function of c, with (a,b) fixed in standardized form."""
if c <= 0:
return -np.inf
x = np.asarray(x, dtype=float)
if np.any((x <= a) | (x > b)):
return -np.inf
logZ = _logZ_truncweibull_min(c, a, b)
n = x.size
return float(n * np.log(c) + (c - 1.0) * np.sum(np.log(x)) - np.sum(x**c) - n * logZ)
# Example: MLE for c when (a,b) are known
true_c, true_a, true_b = 1.8, 0.2, 2.0
data = truncweibull_min.rvs(true_c, true_a, true_b, size=2000, random_state=rng)
def nll(c: float) -> float:
return -loglike_c_given_ab(c, data, true_a, true_b)
res = optimize.minimize_scalar(nll, bounds=(0.1, 10.0), method="bounded")
c_hat = float(res.x)
c_hat, true_c
(1.7778192822675947, 1.8)
7) Sampling & Simulation#
Because we have a closed-form PPF, sampling is straightforward via inverse transform.
From the truncated CDF, one convenient identity is
So the algorithm is:
Sample \(U\sim\mathrm{Unif}(0,1)\).
Compute \(S = (1-U)e^{-a^c} + U e^{-b^c}\).
Return \(X = (-\log S)^{1/c}\).
This uses only log, exp, and power, so it is easy to implement with NumPy.
def truncweibull_min_rvs_numpy(
c: float,
a: float,
b: float,
size: int,
rng: np.random.Generator | None = None,
) -> np.ndarray:
"""NumPy-only sampling for standardized truncweibull_min(c, a, b)."""
if rng is None:
rng = np.random.default_rng()
u = rng.random(size)
return truncweibull_min_ppf(u, c, a, b)
c2, a2, b2 = 1.7, 0.2, 2.0
s_numpy = truncweibull_min_rvs_numpy(c2, a2, b2, size=20000, rng=rng)
s_scipy = truncweibull_min.rvs(c2, a2, b2, size=20000, random_state=rng)
np.mean(s_numpy), np.mean(s_scipy), truncweibull_min.mean(c2, a2, b2)
(0.8850091193506309, 0.8853549491974178, 0.8838055783319058)
8) Visualization#
We’ll visualize:
PDF (analytic) and histogram of Monte Carlo samples
CDF (analytic) and empirical CDF
c3, a3, b3 = 1.7, 0.2, 2.0
samples = truncweibull_min_rvs_numpy(c3, a3, b3, size=5000, rng=rng)
x = np.linspace(a3 + 1e-6, b3, 600)
pdf = truncweibull_min_pdf(x, c3, a3, b3)
cdf = truncweibull_min_cdf(x, c3, a3, b3)
fig = make_subplots(rows=1, cols=2, subplot_titles=("PDF + histogram", "CDF + ECDF"))
# PDF + histogram
fig.add_trace(go.Histogram(x=samples, histnorm="probability density", nbinsx=40, name="samples"), row=1, col=1)
fig.add_trace(go.Scatter(x=x, y=pdf, mode="lines", name="pdf"), row=1, col=1)
# CDF + ECDF
ecdf_x = np.sort(samples)
ecdf_y = np.arange(1, ecdf_x.size + 1) / ecdf_x.size
fig.add_trace(go.Scatter(x=x, y=cdf, mode="lines", name="cdf"), row=1, col=2)
fig.add_trace(go.Scatter(x=ecdf_x, y=ecdf_y, mode="lines", name="ecdf"), row=1, col=2)
fig.update_xaxes(title_text="x", row=1, col=1)
fig.update_xaxes(title_text="x", row=1, col=2)
fig.update_yaxes(title_text="density", row=1, col=1)
fig.update_yaxes(title_text="probability", row=1, col=2)
fig.update_layout(title=f"truncweibull_min(c={c3}, a={a3}, b={b3})", bargap=0.05)
fig.show()
9) SciPy Integration#
scipy.stats.truncweibull_min provides:
pdf,logpdf,cdf,ppf,rvsstatsandmomententropyfit(MLE-based parameter estimation)
A common workflow is to freeze the distribution with parameters and then call methods on the frozen object.
c4, a4, b4 = 1.8, 0.2, 2.0
rv = truncweibull_min(c4, a4, b4) # frozen standardized distribution
x0 = np.array([0.25, 0.5, 1.0, 1.8])
rv.pdf(x0), rv.cdf(x0), rv.ppf([0.1, 0.5, 0.9])
(array([0.597198, 0.847306, 0.72325 , 0.176502]),
array([0.027815, 0.213994, 0.631769, 0.972303]),
array([0.357434, 0.830846, 1.51064 ]))
# Fitting: often (a,b) are known measurement bounds, so we fit c with a,b fixed.
sim = rv.rvs(size=3000, random_state=rng)
c_hat_fit, a_hat_fit, b_hat_fit, loc_hat_fit, scale_hat_fit = truncweibull_min.fit(sim, fa=a4, fb=b4, floc=0, fscale=1)
c_hat_fit, (a_hat_fit, b_hat_fit, loc_hat_fit, scale_hat_fit)
(1.787695312500002, (0.2, 2.0, 0, 1))
10) Statistical Use Cases#
(a) Hypothesis testing#
A simple and interpretable test is whether the data are compatible with an exponential tail on \((a,b]\). Because Weibull with \(c=1\) is exponential, we can test
We can use a likelihood ratio test (LRT) when \((a,b)\) are fixed/known.
(b) Bayesian modeling#
When you want full uncertainty quantification, treat \(c\) as a parameter with a prior (e.g. a Gamma prior) and compute the posterior.
(c) Generative modeling#
Once parameters are estimated, the distribution is a handy bounded generator for synthetic lifetimes or severities (e.g. simulation studies, stress testing, Monte Carlo pipelines).
# (a) Likelihood ratio test for H0: c=1 vs H1: c free (a,b known)
a5, b5 = 0.2, 2.0
x_obs = truncweibull_min.rvs(1.8, a5, b5, size=800, random_state=rng)
c_mle = float(optimize.minimize_scalar(lambda c: -loglike_c_given_ab(c, x_obs, a5, b5), bounds=(0.1, 10.0), method="bounded").x)
ll_alt = loglike_c_given_ab(c_mle, x_obs, a5, b5)
ll_null = loglike_c_given_ab(1.0, x_obs, a5, b5)
lrt = 2.0 * (ll_alt - ll_null)
p_value = float(stats.chi2.sf(lrt, df=1))
dict(c_mle=c_mle, lrt=lrt, p_value=p_value)
{'c_mle': 1.904282724745143,
'lrt': 83.81779762028248,
'p_value': 5.425417704887916e-20}
# (b) Simple Bayesian inference for c with a,b fixed: grid posterior
c_grid = np.linspace(0.2, 5.0, 600)
# Prior: Gamma(shape=2, scale=1) on c (mean=2). Feel free to adjust.
log_prior = stats.gamma(a=2.0, scale=1.0).logpdf(c_grid)
log_like = np.array([loglike_c_given_ab(c, x_obs, a5, b5) for c in c_grid])
log_post_unnorm = log_prior + log_like
log_post = log_post_unnorm - special.logsumexp(log_post_unnorm) - np.log(c_grid[1] - c_grid[0])
post = np.exp(log_post)
c_mean = float(np.trapz(c_grid * post, c_grid))
# 95% credible interval by CDF on the grid
cdf_post = np.cumsum(post) * (c_grid[1] - c_grid[0])
c_lo = float(np.interp(0.025, cdf_post, c_grid))
c_hi = float(np.interp(0.975, cdf_post, c_grid))
fig = go.Figure()
fig.add_trace(go.Scatter(x=c_grid, y=post, mode="lines", name="posterior"))
fig.add_vline(x=c_mean, line_dash="dash", annotation_text=f"mean={c_mean:.2f}")
fig.add_vrect(x0=c_lo, x1=c_hi, opacity=0.15, annotation_text="95% CI")
fig.update_layout(title="Posterior over c (a,b fixed)", xaxis_title="c", yaxis_title="density")
fig.show()
dict(c_mean=c_mean, c_95=(c_lo, c_hi))
{'c_mean': 1.897811196254147, 'c_95': (1.7183616784367284, 2.0639336197638096)}
# (c) Generative modeling: fit c, then generate synthetic samples
c_fit = float(truncweibull_min.fit(x_obs, fa=a5, fb=b5, floc=0, fscale=1)[0])
rv_fit = truncweibull_min(c_fit, a5, b5)
synthetic = rv_fit.rvs(size=5000, random_state=rng)
fig = px.histogram(
x_obs,
nbins=40,
histnorm="probability density",
opacity=0.55,
labels={"value": "x"},
title=f"Observed vs synthetic (fit c={c_fit:.2f}, a={a5}, b={b5})",
)
fig.add_trace(go.Histogram(x=synthetic, nbinsx=40, histnorm="probability density", opacity=0.55, name="synthetic"))
x_line = np.linspace(a5 + 1e-6, b5, 600)
fig.add_trace(go.Scatter(x=x_line, y=rv_fit.pdf(x_line), mode="lines", name="fitted pdf"))
fig.update_layout(barmode="overlay", xaxis_title="x", yaxis_title="density")
fig.show()
11) Pitfalls#
Parameter validity: require
c>0,a>=0,b>a, andscale>0.Standardized bounds: in SciPy,
aandbare standardized; the actual cutoffs areloc + a*scaleandloc + b*scale.Data outside support: any observation \(x\le a\) or \(x>b\) has likelihood 0 under the model.
Numerical stability: the normalizer \(e^{-a^c}-e^{-b^c}\) can suffer cancellation when \(b^c-a^c\) is tiny; prefer log-space or
expm1-based forms.Fitting: if you let
fitestimateaandb, it may push them toward sample min/max; in many applications truncation bounds are known and should be fixed.MGF existence when \(b=\infty\): for infinite upper bound, the MGF can diverge for \(t>0\) depending on
c(like the base Weibull).
12) Summary#
truncweibull_minis a Weibull minimum conditioned to lie in \((a,b]\) (witha,bstandardized in SciPy).The PDF is the Weibull PDF divided by \(Z=e^{-a^c}-e^{-b^c}\); the CDF and PPF follow in closed form.
Raw moments reduce to incomplete gamma differences, enabling mean/variance/skewness/kurtosis calculations.
Sampling is easy via the PPF and can be implemented with NumPy only.
For inference,
scipy.stats.truncweibull_minsupports evaluation, simulation, and MLE-based fitting; consider fixing truncation bounds when they are design constraints.